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SUMMARY 


Three dimensional large-eddy simulations (LES) of single and multiple jet-in- 
cross-flow (JICF) are conducted using the 19-bit Lattice Boltzmann Equation (LBE) 
method coupled with a conventional finite- volume (FV) scheme. In this coupled 
LBE-FV approach, the LBE-LES is employed to simulate the flow inside the jet 
nozzles while the FV-LES is used to simulate the crossflow. The key application 
area is the use of this technique is to study the micro blowing technique (MBT) 
for drag control similar to the recent experiments at NASA/GRC. It is necessary 
to resolve the flow inside the micro-blowing and suction holes with high resolution 
without being restricted by the FV time-step restriction. The coupled LBE-FV-LES 
approach achieves this objectives in a computationally efficient manner. A single 
jet in crossflow case is used for validation purpose and the results are compared 
with experimental data and full LBE-LES simulation. Good agreement with data is 
obtained. Subsequently, MBT over a flat plate with porosity of 25% is simulated using 
9 jets in a compressible cross flow at a Mach number of 0.4. It is shown that MBT 
suppresses the near- wall vortices and reduces the skin friction by up to 50 percent. 
This is in good agreement with experimental data. 
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CHAPTER I 


INTRODUCTION 


In recent years, a novel drag reduction method in turbulent boundary layers has been 
demonstrated using a micro-blowing technique (MBT) at NASA GRC [1, 2, 3, 4, 5, 
6]. Key features of this technique (in contrast to earlier more conventional blowing 
method) are the low effective roughness for the porous skin (achieved due to the 
use of micro holes) and the minimal amount of injection flow rate needed to achieve 
drag control. Results under laboratory and full-scale operating conditions show that 
significant skin friction drag reduction (as much as 50-70 percent) can be achieved 
using MBT. MBT appears to show even higher drag reduction in supersonic boundary 
layers with an added advantage of reduction of noise. 

Interesting observations in subsonic flow are that the drag can be controlled by 
adjusting the injection flow rate and that the maximum drag reduction appears to 
occur within regions of adverse pressure gradient. However, more recent results in 
strongly adverse pressure gradient flow on a strut suggest that micro-blowing can lead 
to increased boundary layer and wake thickness [5] which can result in an increase 
in pressure drag for external flows. Thus, there are still many unresolved questions 
regarding the underlying physics of drag reduction as achieved by the MBT and 
how this blowing process impacts the large-scale flow features. Furthermore, the 
experimental data suggests that the injection system (and the injected airflow) couple 
strongly with the outer (primary) flow especially in adverse pressure gradient flow. 
This makes the optimization of the design of this device problematic, and parametric 
study using primarily an experimental approach is not a cost-effective approach. 

Numerical studies of the MBT in subsonic turbulent boundary layers have also 
been recently reported [7] in which a steady-state 3D code was used. Only a very few 
micro holes and the cross-stream boundary layer were numerically modeled using a 
low Reynolds number turbulence closure. Results showed that for all the simulated 
cases, micro-blowing leads to unsteadiness due to the formation of vortices. However, 
due to the steady-state model employed this feature could not be studied. 

MBT studies using unsteady simulation is problematic since many holes (typically 
hundreds) may have to be simulated simultaneously in order resolve the dynamics. 
Moreover, to properly capture this dynamics, the flow inside the holes has to be com- 
puted so that the flow at the hole exit plane evolves naturally (the flow can be blowing 
or suction depending upon the local pressure gradient). This type of requirement im- 
plies a very high computational overhead due to the increased resolution, and also 
due to the reduced time-step if a conventional FV scheme is employed. Clearly, an 
alternate technique is needed if this type of flow physics has to be simulated. 

In the past, MBT studies in supersonic boundary layers were conducted [8] in 
which Micro-Electro-Mechanical System (MEMS) based micro-scale blowing/suction 
devices were simulated using direct numerical simulation (DNS) in 2D. Significant 
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unsteadiness is shown to be associated with the interaction between the injected fluid 
and the supersonic boundary layer. Vortex shedding and even pairing occur in the 
near field of the injected fluid in the boundary layer. By combining blowing and 
suction it was further shown that not only drag can be reduced but also that it can 
be increased at specific location. 

It is clear that DNS of 3D high Reynolds number turbulent boundary layer with 
proper resolution of the micro holes is beyond the current and perhaps future com- 
putational capability. An alternate method that has the potential for such a study 
is based on large-eddy simulations (LES). In LES, all scales larger than the grid axe 
resolved in a space and time accurate manner and only the scales smaller than the 
grid are modeled using a subgrid model. Although LES looks promising, it also has 
a serious problem in near-wall flows. To properly resolve the small-scale dynamics 
in the log layer, the wall normal resolution has to be close to the DNS requirement. 
This implies that the computational cost will be unacceptable. Clearly, an alternate 
method is required. 

Here, we consider a hybrid approach that has the potential for dealing with the 
flow inside and outside the micro holes. In this approach, we combine a LBE-LES 
method with the conventional FV-LES method to achieve proper resolution in both 
regions. Both 2D and 3D LBE-LES methods have been demonstrated for fuel-air 
mixing and jet in crossflow in the past [9, 10, 11, 12]. In all these earlier reported 
studies, the flow outside and inside the injector or jet inlet are resolved completely 
using LBE. Thus, the flow at the jet exit plane evolves in a natural manner. The 
unique feature of the LBE approach is that it solves the Boltzmann equation (which, 
in the continuum limit recovers the Navier-Stokes equation). Since the Boltzmann 
equation is a single scalar equation, it is computationally very efficient (in fact, orders 
of magnitude faster than conventional FV algorithm). Thus, very high resolution (in 
fact, DNS-like resolution) can be use in the LBE model without a significant increase 
in the computational cost. 

In the present study, the LBE model is employed primarily to resolve the flow 
field inside the injectors while the conventional FV-LES model is used to simulate 
the boundary layer flow. This approach takes the best of both worlds and couples 
them together within a single formulation. The LBE solver is fully coupled to the 
LES solver and interacts across block structured grid domains. Thus, in the injection 
port regions and in the near- wall region (e.g., below y + = 100), a high (DNS-like) 
resolution can be used, and the flow field simulated without requiring any modeling 
using the LBE model, while in regions away from the wall, a conventional FV LES 
code can be employed. 

This report describes the development and the evaluation of the fully coupled 
LBE-LES and FV-LES simulation of turbulent jet in crossflow and multiple micro- 
jets in crossflow. Earlier, the single free jets and jet in crossflow (JICF) were simulated 
using a LBE-LES method [12, 13, 10, 14, 11] and compared to experimental data. 
Very good agreement was obtained. Here, we simulate the same test case using the 
coupled formulation and compare to earlier studies to demonstrate the applicability 
of this coupled method. Subsequently, the coupled solver is used to simulate a matrix 
of 3 x 3 micro holes in a high Mach number turbulent boundary layer. This test case 
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is similar to the experiment at NASA/GRC [4, 6] except that only a 3 x 3 matrix is 
simulated, whereas in the experiment, thousands of micro holes on a large flat plate 
were employed. 

Here, we demonstrate this coupled approach and show the feasibility to extend 
this technique to handle many holes simultaneously. We also believe that these sort of 
sub-set simulations can be used to investigate similarity features and also to develop 
new scaling laws. However, due to resource and time constraints, some of the planned 
tasks had to be truncated. Nevertheless, as this report demonstrates, we now have a 
hybrid LBE-FV capability to carry out LES of this complex problem. 

We have reported on our past progress elsewhere in some detail [12, 13, 10, 14]). 
Therefore, the details provided in these cited papers/reports are not included in this 
report, for brevity. Only new issues addressed subsequently are discussed in this 
paper. We describe the issues addressed in developing the hybrid solver and some 
of the validation studies. Finally, we show how this code can be used to simulate 
an array of micro-hole blowing into a compressible crossflow. We believe we have 
achieved all the key objectives of this research but full exploitation of this capability 
remains to be fully explored. 
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CHAPTER II 


SIMULATION MODEL 


The simulation model for the finite-volume and lattice approaches are summarized 
here for completeness. More details are given in the cited references. 

2.1 FV-LES Formulation 

Our interest here is the application of MBT to reduce skin friction drag on air- 
craft wings under realistic flight conditions in the Mach number range of 0.4 - 1.8. 
Therefore, we employ a fully compressible Navier-Stokes solver that has been used 
extensively in the past and is fully validated in non-reacting [15] and reacting flows 
[16, 17, 18]. The solver is nominally second-order accurate in both space and time; 
however, it also has the option for fourth order accuracy in space. 

The LES equations are obtained by the application of a spatial filter to the gov- 
erning equations of motion. For the compressible equations, a Favre spatial filter is 
employed to separate the resolved and unresolved motion [19]. For a finite-volume 
scheme, a low-pass, top-hat filter of the local grid size (A) is appropriate. The filtering 
operation results in the compressible LES equations. 

Applying the filtering to the conservation equations results in the following LES 
equations: 


dp dpui 
dt dxi 


(2.1) 


+ + p6ij ~ fij + = 0 ( 2 - 2 ) 

_ Ujfji + H- 9S + cq 9S ] = 0 (2.3) 

In the above equations, all terms with superscript “sgs” denote subgrid quanti- 
ties that require closure. Here, p and iq are respectively, the mixture density and 
the velocity. Also, is the filtered viscous tensor (defined in terms of the filtered 
quantities) and q t is the filtered heat flux vector given by 

In the above relations, p and 8; are respectively, the filtered molecular viscosity and 
the thermal conductivity. These transport properties are determined as a function of 
the filtered quantities. 

The pressure is determined from the filtered equation of state, p = pRT and the 
filtered total energy per unit volume is pE = pe + \puiUi + pk S9S . Here, R is gas 
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constant and the subgrid kinetic energy is defined as, k sgs = 
Finally, e is filtered internal energy for calorically perfect gas. 

(l/2)[u fc itfc - u k u k }. 

The subgrid terms that require closure are: 


P {UiUj Ui Uj ) 

(2.5) 

H° 9S = p (Eui - Etii) + (pul - put) 

(2.6) 

(Tj = UjTij — UjTij 

(2.7) 


2.1.1 Momentum transport closure 

The subgrid stress tensor t*?* is closed using an subgrid eddy viscosity and a gradient 
diffusion model at the grid cut-off scale: 

r‘f s = -2 putlSij - IslkSy] + (2.8) 


Here, the resolved rate of strain is given as = (1/2 )(dui/dxj + dilj/dxi ) and 
the subgrid eddy viscosity is v t . Since large-scale motion is resolved in LES, the 
associated counter-gradient effects at the large scales are also resolved (even though 
a gradient closure is employed for r- 9S at the grid cut-off scale). 

To obtain u t and k sgs , a non-equilibrium model [20, 21] for the subgrid kinetic 
energy, k sgs is used, and is given by : 


♦£<**->--** I 


dt 


r Jk a 9 S f /2 

Cip ~A~ 


+ 


d {_v t dk aga \ 


dxi \ C t dxi ) 


(2.9) 


The subgrid eddy viscosity is modeled as: v t = C u A\/k S9S , where A = (AxAyAz) 1 ^ 3 
is based on local grid size (Ax, Ay, A z). There are three coefficients in the above clo- 
sure coefficients, C u , C t , and C t . The nominal “constant” values for these coefficients 
are [22, 23] 0.067, 0.916 and 1.0, respectively. 

However, these coefficients can be obtained by using a localized dynamic procedure 
for the subgrid kinetic energy model[16, 24] (denoted as LDKM, hereafter). This 
dynamic procedure uses the experimental observation in high Re turbulent jet [25] 
that the subgrid stress t*? s at the grid filter level A and the Leonard’s stress J Ly(= 

[{pUiUj) — ± p ^ -~ ) at the test filter level A(= 2A) are self-similar. Here (and 

henceforth), < / > and / both indicate test filtering. Since Ly can be explicitly 
computed at the test filter level, a simple scale-similar model of the form t-? s = C^L y , 
where Cl is an adjustable constant, was proposed earlier [25] but was found to lack 
proper dissipation. 

In the LDKM model, the above observation is extended and it is assumed that Ly 
and the subgrid stress rf/ s at the test filter level are also similar (i.e., Ly = CiT^f 3 )- 
Using this, r//' s is modelled using the same form as for (Eqn. 2.8), except that 
all variables are defined at the test_filter level. We define the subgrid kinetic energy 

at the test filter level as k test = 5 [^ — ^h] and obtain a relation 

* p p 

Lij = CluT = -2^av^A((4) - ^(Su)<5y) + 2 -ClP k tes A 3 (2.10) 
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In the above equation, we assume Ci = 1 and so, the only unknown is C v . This 
equation is, thus, an explicit model representation for the constant C„ in terms of 
quantities resolved at the test filter level. This system of equations represents five 
independent equations for one unknown coefficient (and hence, is an over-determined 
system). The value of C„ is determined in an approximated manner by applying the 
least-square method [26] . Thus, 


In the above expression 


L'nMn 

2 MijMij 


Lij — Lij ^ P ktest^ij 


( 2 . 11 ) 


( 2 . 12 ) 


^ , 1 ~ 

M ij p V ktest^{{Sij) g (Skk)fiij^ 


(2.13) 


A similar approach is used to obtain the dissipation coefficient C e such that: 


^ _ A + 


ilj 


ohtj 

dxi 


(2.14) 


where p is the molecular viscosity and p t (— v t * p) is eddy viscosity at the grid filter 

level. The tensor T i} is defined as and T t j indicates tensor at the 

test-filter level. 

More details are given elsewhere[16, 24]. There are a few noteworthy points to 
highlight in this closure: (a) the LDKM approach does not employ the Germano’s 
identity [27], (b) the self-similar approach implies that both A and A must lie in 
the inertial range, and this provides a (albeit) rough estimate for the minimum grid 
resolution that can be used for a given Re, (c) the denominator in both Eqns. (2.11) 
and (2.14) are well-defined quantities at the test filter level, and can be directly 
computed, (d) the evaluation of the coefficients can be carried out locally (i.e., at all 
grid points) in space without encountering any instability, (e) the LDKM approach 
satisfies all the realizability conditions [28] in the majority of the grid points even 
in complex swirling reacting flows, and (f) the dynamic evaluation can be used near 
walls without any change [29, 30]. Finally, the computational overhead of the LDKM 
is not very significant since only one additional equation has to be solved. 


2.1.2 Energy and scalar transport closure 


In addition to r?f a , other subgrid terms that appear in the LES filtered energy and 
species equations have to be closed. The subgrid total enthalpy flux, H- 9S is also 
modeled using the eddy viscosity and a gradient assumption as: 


Hl gs = 


P Pr t dxi 


(2.15) 
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Here, H is the filtered total enthalpy and Pr t is a turbulent Prandtl number that 
can also be computed using a dynamic procedure but is currently assumed to be 
unity. The total enthalpy term H is evaluated as sum of specific_enthalpy of mixture, 
specific kinetic energy, and specific sub-grid scale energy: H = h + ^ 4- k sgs . 

The other unclosed terms such as a° 9S are often neglected in conventional LES 
approaches [31] but there is no clear justification. This issue is not addressed in this 
study. 

2.2 LBE-LES using Single Relaxation Time 

The LBE method originates from a Boolean fluid model known as the lattice gas au- 
tomata (LG A) which simulates viscous fluid flow by tracing the fluid motion through 
advection of fluid particles and particle collision on a regular lattice. LBE is an 
improvement over LGA in which the Boolean fluid model is replaced by a single con- 
tinuous particle distribution, which is analogous to the particle distribution function 
in kinetic theory. This replacement eliminates the intrinsic noise inherent in LGA 
schemes and overcomes the shortcomings of a limited transport coefficient. The in- 
troduction of the BGK single relaxation time (SRT) model for the collision operator 
further simplifies the algorithm and eliminates the lack of Galilean invariance and 
the dependence of pressure on velocity [32, 33]. This model assumes that the parti- 
cle distribution function relaxes to its equilibrium state at a constant rate, and the 
collision operator is similar to the classical BGK Boltzmann operator [34]. 

Whereas conventional Navier-Stokes schemes solve the macroscopic properties of 
the fluid explicitly, LBE method solves the Boltzmann equation by tracking the evo- 
lution of the microscopic particle distribution of the fluid in phase space (velocity 
space, physical space and time). Consequently, the conserved variables of the fluid 
(density and momentum) are obtained indirectly by local integration of the particle 
distribution (over the velocity space). The incompressible Navier-Stokes is recovered 
in the nearly incompressible limit of LBE using the Chapman- Enskog expansion. For 
the present MBT application this is an acceptable model since we are interest in using 
the LBE approach only within the jet injectors where the flow is at a very low speed. 

Solving the lattice Boltzmann equation instead of the Navier-Stokes equation pro- 
vides three distinct advantages. First, due to the kinetic nature of the LBE method, 
the convection operator is linear. Simple convection in conjunction with a collision 
process allows the recovery of the nonlinear macroscopic advection through multi- 
scale expansions. Second, because the macroscopic properties of the flow field is not 
solved directly, LBE method avoids solving the Poisson equation, which proves to 
be numerically difficult in most finite difference methods. Third, the macroscopic 
properties are obtained from the microscopic particle distributions through simple 
arithmetic integration. More details are given in a recent review [35]. 

LBE method consists of two primary steps. The particles first stream to its next 
nearest neighbor in the direction of its prescribed velocity. Subsequently, particles of 
different velocities arriving at the same node interacts with each other by relaxing to 
its local equilibrium values which are formulated specifically to recover the low Mach 
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number limit of the Navier-Stokes equation. The evolution of the non-dimensional 
distribution function f a is thus governed by: 


/ a (x + e a <5,t + <5) -/ a (x,t) = -[/^(x,t) - / a (x,f)], 

T 

a = 0,1, •••,18 (2.16) 

where r is the relaxation time, is the equilibrium distribution function and e a is 
the particle speed in a direction. The characteristic speed is thus c = e a 5/5 = |e a |. 
Rest particles of type 0 with eo = 0 are also allowed. Note that the time step and the 
lattice spacing each have equal spacing of unity. Thus, <5 = 1. More recently, studies 
have demonstrated that multiple relaxation time (MRT) models can be developed 
when one needs to simulate scalar mixing at different Schmidt number [36]. The 
MRT approach has also been shown to be more stable in complex flows and in high 
Reynolds number flows. Therefore, we have also reviewed its performance for the 
present study. 

In principle, there are an infinite number of possible velocity directions in the 3D 
velocity space. Discretizing these infinite number of velocity directions into a fixed 
set of velocity directions inevitably introduces discretization errors to the solution. 
As a general rule, the accuracy of the model to simulate Navier-Stokes flow comes at 
the expense of increasing computational cost resulting from the number of discrete 
velocities used in the model. Frisch et al. [37] have shown that the Navier-Stokes 
equation cannot be recovered unless sufficient discrete velocities is used to ensure 
lattice symmetry. 

There are various 3D cubic lattice models developed, most notably the 15-bit 
(D315), 19-bit (D3Q19), and 27-bit (D3Q27) model [38]. Here, using common nota- 
tions in scientific literatures, D is the number of dimensions and Q is the number of 
discrete velocities. In previous numerical simulations of a square duct, a lid-driven 
cavity and a circular pipe [39], no significant improvement in accuracy is observed 
when the D3Q27 model was used over the D3Q19 model, and thus, the D3Q19 model 
is assumed to be sufficiently accurate for the current purpose. 

The 19-bit velocity field (Fig. 2.1) is: 


e a = < 


f (0,0,0) 

for a = 0, rest particle, 
((±l,0,0),(0,±l,0),(0,0,±l))c 
for a = 1, 2, • • • , 6, class I links, 

((±1, ±1, 0), (0, ±1, ±1), (±1, 0, ±1 ))V2c 
^ for a = 7, 8, • • • , 18, class II links. 

Here, / a 9 is given by the following form: 


fa — w aP[ 1 + 


3(e a • u) 9 (e a • u) 2 3u 2 . 


+ ^ 


2 c 2 


where 


j a = 0 

u>a= { i a = 1,2, ••• ,6 

55 a = 7,8, ,18. 


(2.17) 


(2.18) 
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The macroscopic properties of the flow field can be obtained by integrating the 
distribution functions over the velocity space: 

P = (2-19) 

a 


pu = J% Q / Q , a = 0,1,..., 18 (2.20) 

a 

where p is the density and u is the velocity. All the macroscopic properties are 
obtained as function of space and time from this integration. 

The Navier-Stokes mass and momentum equations are obtained by using the BGK 
single relaxation time model [34], and by employing Chapman- Enskog expansion are: 

dp dpui 
dt + dxi 



d(p Uj) dpUjUj = d{c 2 s p) fflvpSjj . . 

dt dxj dxi dxj 

Here, repeated indices indicate summation and S l3 = |( dui/dxj + duj/dxi) is the 
strain-rate tensor. The non-dimensional pressure is given by the constant temperature 
ideal gas equation of state p — c^p where c s is the speed of sound with (c s = c/V 3), 
and v = [(2 r — l)/6] is the kinematic viscosity. 

For LES application we consider the filtered distribution function such that the 
LBE-LES form of the governing equation is 

/a(x + e a 6, t + 6 ) - ja(x, t) = — [/a 9 (x, t) ~ /«(x, t)] 

T sgs 


a = 0,1, ■■■,18 (2.23) 

where the distribution function f a represents only those of the resolved scales. The 
effect of the unresolved scale motion is modeled through an effective collision term, 
which in the BGK approximation is included as an effective relaxation time T sgs . The 
form of the subgrid correction is not fully explored at present. 

In the earlier studies [12] the effective relaxation time was obtained using an eddy 
viscosity model (based on the Smagorinsky’s eddy viscosity model) such that: 

•> + (2.24) 

6 

with the eddy viscosity u T determined using: 

v T = GA 2 5 (2.25) 


where C v is the Smagorinsky constant, A = (A x A y A 2 )5 is the associated length scale 
and S = | S tJ Sij | . Here, Si 3 = | (<9n, jdx 3 + du 3 / dx x ) is the resolved-scale rate-of-strain 
tensor is the characteristic filtered rate of strain tensor. 
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The Smagorinsky constant C v is determined using the localized dynamic model 
(LDM) [40, 24] described earlier but now adapted for the Smagorinsky’s model (this 
is done primarily for convenience at this time). Thus, C v is obtained using 


2 MijMij 

where Lk = L tj - lfjk test S i:j and M tj = k test A((Sij) - l(S~ kk )5 ij ). Here, L {j = 
J)((uiUj) — ( Ui)(uj )) is the Leonard stress tensor and k test = \{{uiUj) — ( vk)(uj )) = 
\L kk fp is the resolved kinetic energy at the test-filter level. 

This algebraic model has some inherent limitations in conventional LES and these 
limitations remain in effect in the LBE-LES approach. Thus, for high Re flows, the 
lattice resolution has to be quite high since the algebraic model requires that nearly 
all of the inertial range is resolved. The alternate one-equation model [24] has the 
ability to deal with high-Re flows on relatively coarse grid and therefore, is should be 
implemented for future studies. 

Finally, it is worth noting that this type of closure also has some fundamental 
limitation in the LBE context since it is based on a model in the physical space, 
whereas the LBE simulation is in the phase space. A proper subgrid representa- 
tion for the LBE-LES will require a closure model in the phase space to model the 
subgrid component of the collision integral. Such an approach still remains to be 
demonstrated. 



2.3 LBE-LES using Multiple Time Relaxation 

The simplest collision model is the linearized BGK model using single relaxation 
time (SRT) for all velocity bits described in the previous section. Since this model 
has some restrictions as well as stability issues near the sharp edges, Lallemand et. al 
[41] developed a generalized LBE model using multiple relaxation time (MRT). MRT 
model is constructed in moment space and it characterizes the instability, dissipation 
and dispersion of the LBE model. In addition, it provides the maximum number 
of adjustable parameters which can be used to optimize its properties. This section 
covers the formulation of the MRT. 

The evolution of the non-dimensional distribution function f a based on MRT 
model is governed by: 

[/a( x + e cA) t + <5i)](jV+l)xl — [/q( x , t)](W+ l)xl = 

- [^^(iv+^xCiv+qit/QCx^jl^+qxi - [fl q (x,t)\( N+ i) xl } 

where f* 9 is the equilibrium distribution function for each velocity link of a, which 
goes from 0 to N. The velocity direction for each link is indicated by e a with the 
magnitude of |e Q |, which is the normalized characteristic speed of particles. As in 
the SRT model, the speed for rest particles is equal to zero (e 0 = 0) and for class I 
particles is equal to ratio of grid size to time step (|e Q | = S x /S t = 1). Note that, all 
time and lengths in LBE formulation are normalized by the time step and the grid 
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size respectively, therefore, the normalized time step (S t ) and grid size (S x ) themselves 
will be unity. S is the collision matrix and is designed based on its eigenvector (M) 
and eigenvalue (S) matrices. The eigenvector matrix M is called moment matrix and 
it consists of 0 orthogonal basis vectors constructed by Gram-Schmidt orthogonality 
of the velocity links. In this study, D3Q19 link configuration is used. For this link 


configuration , 

the 

moment 

matrix is 

equal to : 
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(2.27) 

Note that this eigenvector matrix is slightly different than what is given in Ref. 
[42] since the velocity links are numbered differently. Another important parameter, 
that can be found by this matrix, is the moments (m) of the distribution function 
using m = Mf . These moments are defined as: 

m = [/£>, e, e, jx , Qx i jy > Qyi jz > Qz ; 3 0" XXi ^^xx > &ww > ^wui > @ xy t @yzi & zx i ^2-xi rn 2 ] 

where p, e, e, j, q, a, and m are respectively, density, kinetic energy, kinetic energy 
square, momentum, energy flux, traceless viscous stress tensor (with a ww = cr yy — (T zz ), 
and antisymmetric third-order moment. The moments tt ix and n ww are proportional 
to a xx and a ww with the factor of w xx which is an additional free parameter for D3Q19 
model. This parameter is equal to zero for optimized stability and is equal to -0.5 to 
recover the LBGK model. 

The eigenvalue matrix is a diagonal matrix as following: 

S = diag(Ao, Ai, A2, • • • , Ais) 

The diagonal elements (A a ) of this matrix are equal to the collision frequency of 
the correspondent moment. Such frequencies for conserved moments (hydrodynamic 
variables) of density and momentum is equal to zero [42], so A 0 = A3 = A5 = A7 = 0. 
The collision frequency for viscous stresses are found based on the Chapmann-Enskog 
analysis, which gives: A 9 = An = A J3 = X u = A15 = 1/r = 2/(6i/ + 1). The rest of 
the frequencies can be found based on linear analysis to achieve optimized stability 
for the model, following [42], They are equal to: 

Ai = 1.19, A 2 = A 10 = A 12 = 1.4, 

A 4 = Ag = As = 1.2, 

Al 6 := A 17 = Ais = 1.98 
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The collision term of the generalized LBE (Eqn. 2.3) can be written in terms of 
the moments of distribution function as, 

[/q(x + e a 5, t + £)](jv+i)xi — [/ Q (x, £)](jv+ i)xi = 

— [-^a^ 1 ](W+l)x(Ar+l)[«S , Q/3](Ar+l)x(W+l){[uT. Q (x, t)](w+l)xl — [’’^(x, ^)](AT+l)xl} 

The equilibrium value of conserved moments have the same value as of their pre- 
collision moments, whereas non conserved moments are found as a function of the 
conserved moments. Scalar non-conserved moments are functions of (p, \j\ 2 ) and 
vector moments are functions of (pj, j\j\ 2 ). These function are given in Ref. [42], 

All other issues such as the 19- velocity model and the subgrid closure issues remain 
unchanged. 

2.4 Boundary Conditions 

Many types of boundary conditions have to be implemented for both the FV-LES 
and the LBE-LES algorithms. 

2.4.1 Boundary Conditions for FV-LES 

The compressible FV-LES solver requires characteristic inflow and outflow boundary 
conditions [43] in order to maintain physical accuracy of incoming and outgoing pres- 
sure waves. Details are given elsewhere [43] and therefore, not repeated here. Some 
details related to these BCs for the hybrid solver are discussed later. On physical 
walls are no-slip adiabatic conditions are employed for all these simulations. 

2.4.2 Boundary Conditions for LBE-LES 

For the LBE solver, proper implementation of boundary conditions is very impor- 
tant, especially for complex geometries. On no-slip walls we employ the bounce back 
scheme, i.e., the particles arriving at the stationary wall is reflected back. Although 
easy to implement, the exact location of the no-slip wall is important. The bounce- 
back boundary condition exhibits second order accuracy only when the wall is place 
exactly halfway between the boundary node and the first fluid node [11], Addi- 
tionally, the implementation of no-slip conditions for complex wall geometries (e.g., 
convex and/or concave walls) are also complex. Analysis of these conditions were 
reported earlier but are included here for completeness. 

The convention of the link vectors used in the examples is illustrated in Fig. 2.1 
and their surface orientations in Fig. 2.2. 

For the planar node, if a normal velocity w bc is prescribed, the tangential momen- 
tums of the parallel links and the density are: u tan = (/j + / 7 + / 9 ) — (/ 2 + / 8 + / 10 ) 
Vtan — (/3 + h + /lO) — (fi + fs + / 9) 

P = ( 2 (/e + /12 + /13 + /i 6 + /17 + /o + /l + fi + Jz + f\ + fi + / 8 + fg + /10)/ (1 — Wbc)- 
The unknown populations are therefore: 
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fa = fa + pWbc 
fn /12 0.5ttt on 

fl4 /l3 4” 0.5Ui O 7i 

fib fib 0. OVtari 

fl8 ~ fl7 + 0.5v tan . 

For nodes on concave edge, the tangential momentum is v tan = fa — f 4 and the 
unknowns are: 

fa=fa 
fa = fa 

fn — /12 

fa — fa ~ 0.25u tan 
fa — fw + 0.25u tan 
fab = fib ~ 0.25 Vtan 
fl8 — fl7 + 0.25l ’tan- 

The buried links (links that do not participate in streaming but in collision) are 
averaged and redistributed. Thus, / 13 = / 14 = 0.5(/i3 + fu) 

For the concave corner, all the class I links exist in known/unknown pairs and 
thus: fa = fa, fa = fa, and fa = fa and 

fa — fa 

fn = /12 
fib = fab 

fa = fao = 0.5(/g + fao ) 

/l3 = /l4 = 0.5(/i3 + fu) 
fl7 — /l8 = 0.5(/i7 + /is). 


Maier has pointed out that convex edges and convex corners cannot satisfy the 
no-slip condition with this scheme because there are insufficient unknown populations 
(at least one class I link) at these boundaries [44]. New treatments are devised for 
these boundaries and the no-slip condition is forced through the modification of the 
known populations as well as the unknowns. In short, the unknown populations are 
computed using bounce back, and the known non-parallel populations of opposite 
links averaged and redistributed. Subsequently, the momentum of the parallel links 
is redistributed to the non-parallel pairs. With this method, the total population of 
the known links is preserved. For the geometry in Fig. 2.2: fai = / 12 and 
Vtan f 3 f 4 

fa = fa = 0.5(/i + fa) 
fa = fa — 0 . 5(/5 + fa) 

fl3 = fl4 — 0 - 5 (/ i 3 + fu) 

fa = 0.5(/ 7 + fa) - 0.125w tan 
fa = 0.5(/ 7 + /g) + 0.125u ton 
fa = 0.5(/q + fao) + 0.125i> tan 
fao = 0.5(fa + fao) — 0.125 Vtan 
fib — 0.5(/i5 + fao) — 0. 125Ufan 
fib = 0.5(/i5 + fao) + 0.125u tan 
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fn — 0.5(/i7 + /is) — 0.125vt an 
/is = 0.5(/i7 + /ig) + 0.125u tan . 


For the convex corner, there are no unknowns and therefore all the opposite link pairs 
are averaged and redistributed. 

For the symmetric boundary condition, the unknown populations beyond the 
plane of symmetry are needed during streaming. Thus, the outgoing populations 
are reflected back into the domain exactly the opposite way they leave. For the ge- 
ometry in Fig. 2.2, using k to denote the known outgoing populations at the i plane 
and uk to denote the incoming unknowns at i + 1 , the conditions are: 

ruk } i+l rkj 

J 2 “ J 1 

rufc,i+l rk,i 

J8 — J 9 

ruk,i -\- 1 rk,i 

J 10 ~ J 7 

ruk,i+ 1 rk,i 

J 12 — / 13 

ruk } i+ 1 rkyi 

J 14 “ /ll * 


2.5 LBE-FV coupling 

The two methods, LBE-LES and FV-LES have to be coupled together with full 
two-way interactions in order to solve the problem. Here, we achieve this coupling by 
combining the boundary conditions for FV and LBE models described in the previous 
sections. Two regions need to be properly addressed: (a) a situation when the flow is 
going from the FV-LES domain to the LBE-LES domain and (b) the reverse situation 
of flow from LBE-LES domain into FV-LES domain. In both situations inflow and 
outflow conditions have to be prescribed. In the present study, we are primarily 
interested in micro-blowing and therefore, this allows us to use some simplifications 
(these simplifications can be relaxed for generalize application but not addressed at 
present). Since the diameter of the micro hole is very small (0.5 mm) and the flow 
rate is also very small, we assume that the flow out of these holes is controlled by 
the pressure in the outer domain. Thus, when the pressure is low above the holes 
(relative to the pressure inside the injectors, which is determined primarily by the 
stagnation pressure in the injection system) then there will be outflow. On the other 
hand, if the pressure is higher than the exit pressure of the injector then outflow will 
stop and there will be no net flux. For a general purpose coupling we should allow 
for suction (or inflow into the injectors) to occur in the latter case. Although we 
have looked at this issue, we have not completed this study so far. However, for the 
current application this approach is considered reasonable as a first attempt to couple 
the two methods. 

To provide characteristic “inflow” boundary conditions for the FV-LES (which is a 
compressible solver) one variable must be specified from inside the FV domain while 
four can be prescribed (and one determined by equation of state). In the present 
implementation, temperature at the “inflow” plane is prescribed from inside the FV 
domain (by simple extrapolation) while density and the three velocity components 
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are obtained from the LBE-LES solution from Eqs. 2.19 and 2.20. Finally, pressure 
is computed using the equation of state. This approach is slightly different from the 
approach developed for compressible inflow, as described by Poinsot and Lele [43], 
where velocities and temperature are prescribed and density is computed from the 
backward characteristic. It is possible to update temperature using an appropriate 
characteristic form but has not been addressed here since for all the studies so far, 
the temperature was nearly a constant and simple extrapolation was sufficient, as 
demonstrated later. We plan to revisit this issue in the near future, especially when 
compressibility effects have to be addressed. 

To provide the “outflow” boundary conditions for the LBE-LES from one needs 
to prescribe all the distribution functions f a . Some of these can be determined from 
the FV-LES input and some can be determined from inside the LBE-LES domain 
(see Fig. 2.1 for the lattice directions). In the present study, the pressure at the exit 
plane of the LBE domain is obtained from the FV domain (assuming zero normal 
pressure gradient). Then, the density is computed using p = c 2 s p. The three velocity 
components at the exit plane are given from the LBE domain and then, we need to 
solve six equations in Eq.2.19 and Eq.2.20 along with the symmetry conditions. 

The above noted coupling is spatial in nature but since temporal accuracy is 
needed, this coupling occurs at the proper instant in time. Here, the time scale for 
LBE-LES is obtained by using the characteristic scale of the injector diameter and c. 
This time scale is compared to the characteristic time scale for the FV-LES evolution 
(which is the classical CFL time-step). Both solvers evolve at their individual time 
steps and couple at the appropriate epoch. In the present study, the LBE-LES time 
step is smaller than the FV-LES and thus, the LBE field can evolve longer than the 
FV-LES before they need to be coupled. 


16 



CHAPTER III 


NUMERICAL IMPLEMENTATION 


Most of the details of the numerical implementation has been reported before and 
only some key comments are noted here for completeness. 


3.1 Grid Stretching and Interpolation 

Grid stretching can be carried out in both FV-LES and LBE-LES. For the FV- 
LES, the solver employs a generalized body conforming grid and so grid stretching 
is easy to implement. Typically, stretching less than 5 percent is used in regions of 
interest to maintain second order accuracy. Grid stretching in LBE-LES can also 
be chosen to ensure second order accuracy [11] using the interpolation supplemented 
Lattice Boltzmann (ISLBE). The ISLBE is implemented using the Lagrangian upwind 
quadratic interpolation method. 

For most the reported multi-hole microblowing cases we employed a uniform grid 
in the LBE domain but employed a stretched grid in the FV-LES domain. At present, 
the grid is continuous between the two domains but there is really no need for this 
in this coupled approach. However, to use different grids between the two domains 
interpolation will have to be included. This issue will be addressed in the near future. 


3.2 Grid Resolution Issues 

Grid quality is very important for LES. In general, the numerical dissipation of the 
scheme has to be smaller than the dissipation provided by the subgrid model and 
this implies that both the grid and the numerical accuracy (inherent dissipation) are 
important. Since the current subgrid model is based on the algebraic Smagorinsky 
model, the overall resolution needs to be fine enough to resolve nearly all of the inertial 
range. Some grid resolution effects were addressed in this study using the single jet 
in cross-flow and are reported in the next chapter. 


3.3 Parallel Performance of the FV-LES 

The FV-LES is highly optimized in parallel using the message passing interface (MPI) 
and its performance has be well documented in many papers [45]. This solver runs 
on all parallel systems and maintains high scalability. 
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3.4 Parallel Performance of the LBE-LES 

The LBE-LES solver is also implemented in parallel MPI. Optimization requires 
proper programming of all the boundary conditions and minimizing message pass- 
ing overhead. The computational efficiency of the solver is considerable. For a single 
jet in cross-flow a maximum resolution of 3.25 million grid points have been used for 
the simulations. On Compaq SC45, the LBE-LES solver (with the localized dynamic 
option) costs approximately 400 single-processor hours for 5 flow-through-times and 
with 2.3 GB memory usage. In contrast, a conventional finite-volume solver using 
a similar grid would require substantially more (by 1-2 orders of magnitude) CPU 
hours to complete a similar simulation. 


3.5 Parallel Performance of the Hybrid Solver 

Both the FV-LES and LBE-LES solvers are coupled together domain decomposition 
to achieve proper load balancing. In order to couple the two parallel codes, different 
MPI groups are defined. Thus, the two codes are parallelized and implemented inde- 
pendently and they are coupled through a driver code using message passing. This 
approach allows independent control of both solvers and makes it very easy to use 
LBE-LES in many injector holes without requiring major code revision. Therefore, 
extension of this approach to many hundreds of injectors will not be a major issue 
(other than the computational cost). 
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CHAPTER IV 


RESULTS AND DISCUSSION 


Here, we report on some of the earlier results and new results since the last report. 

4.1 Jet-in- Crossflow (JICF) 

Jet in cross flow (JICF) is a canonical flow that models fuel jet injection in realistic 
devices. In addition, JICF occurs in other real applications, such as, vertical and 
short takeoff/landing (V/STOL) aircraft, internal cooling of turbine blades, dilution 
air jets in combustion chamber of gas turbine engines, etc. In the present context, 
JICF is a generic subset for the MBT approach. Many complex features are observed 
in JICF. In particular, three major structures have been identified in JICF:(l)the 
horse-shoe vortices that form ahead of the jet due to separation of the wall boundary 
layer [46], (2) the vortices that form along the jet shear-layer, and (3) the counter- 
rotating vortex pair (CRVP) that forms in the far wake aligned in the direction of 
the jet [47]. The interaction between these various vortical structures is very complex 
and still not very well understood. For example, it is generally accepted that the 
source of CRVP is the vortices in the jet shear layer but the exact transition from 
one to the other is still under investigation. 

There are several parameters that impact the JICF structures and have been 
studied in the past. Jet to cross flow momentum ratio R and Reynolds number are 
the key parameters [48]. For low-speed single species flows, the velocity ratio can 
be used instead of the momentum ratio. Additional parameters include the effect 
of injection angle, such as the effect of pitched and yeawed injection [49], and the 
effect of the hole exit geometry [50]. As noted by Hwang [6], in the MBT case, 
the micro holes are drilled using laser and all these holes have irregular shapes. We 
limit ourselves to square holes in this study for three primary reasons: (a) the LBE 
approach as implemented here cannot handle curved surfaces, (b) the approximation 
of the irregular micro-holes using a rectangular hole is considered reasonable [6] and 
(c) there is a reasonable experimental data base for a single rectangular jet in cross 
flow [51] with which comparison can be carried out for validation. 

Note that most recently, we have extended the LBE model to deal with curved 
surfaces but this capability has not be fully implemented here. 

4.1.1 Grid Resolution and Conditions 

The experiment of Ajersch et al. [51] is chosen as a benchmark case for validation and 
to investigate in some detail. The dimensions of the computational domain are shown 
in Fig. 4.1. The simulation is carried out at Reynolds number of 4700, based on the 
jet velocity and the nozzle width d. Two cases with jet-cross-flow velocity ratio of 
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Outflow 




Figure 4.1: Geometry and computational domain for the jet in crossflow. 


Grid 

resolution 

time step 

5 

77 

max(^) 

Dense 

50 x 50 x 100 

2.3E-06 

10 

8.0 

Coarse 

26 x 26 x 60 

4.6E-06 

22 

3.0 


Table 4.1: Coarse and dense grid distribution. A is the grid size and 77 is the 
Kolmogorov scale. 


0.5 and 1 are simulated. The cross flow velocity profile is initialized with a boundary 
layer thickness of 2d. The computational domain is resolved using 200 x 150 x 100 
for the crossflow domain and 50 x 50 x 100 for jet section. Thus, a total of 3.25 x 10 6 
grid points is used to discretize the entire domain. Periodic conditions are used in 
the span-wise direction to simulate a single square jet out of a row of six used in 
the experiment. On the top surface slip condition is used and at the exit outflow 
conditions are prescribed. The incoming jet pipe velocity profile with constant value 
is prescribed in the pipe a distance of lOd below the flat plate allowing the flow to 
develop naturally as the jet merges into the crossflow. 

Results are compared with the experiment of Ajersch and numerical results ob- 
tained by Hoda et al [52]. Standard k — e [53], a high- Re, k — e model of Launder- 
Tselepidakis [54], and two-layer turbulence model of Chen [55] have been used in Hoda 
et al studies and therefore, offers an opportunity to compare the current LBE-LES 
results against these different finite- volume closure methods. 

4.1.2 Effect of Grid and Subgrid Model 

Results for a coarse and a dense grid for R = 1 are compared. Table 4.1 gives the 
relevant parameters. 
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(a) Dense grid (y/d = -0.5) (b) Dense grid (z/d = 1.5) 



(c) Coarse grid (y/d = -0.5) (d) Coarse grid (z/d = 1.5) 


Figure 4.2: The ratio of eddy viscosity to physical viscosity at y/d = —0.5 and 
z/d = 1.5 for case R = 1 for (a,b) dense grid and (c,d) coarse grid. 

Figures 4.2 (a)-(d) compare the ratio of eddy and physical viscosity for the coarse 
and dense grid cases. Parts (a) and (b) show the contours at the cross sections of 
y/d = —0.5 and z/d — 1.5, respectively, for the dense grid. Coarse grid viscosity 
ratio contours are shown in parts (c) and (d) for these same locations. The maximum 
viscosity ratio for the coarse grid is about 8, whereas for dense grid is about 3, which 
is qualitatively as expected since eddy viscosity should increase with decrease in grid 
resolution. Furthermore, the eddy viscosity for coarse grid is seen in more regions 
than dense grid, which supports the validity of the localized dynamic model for the 
subgrid model. 

4.1.3 Flow features in JICF 

It is well established in the literature [32, 50] that the jet in crossflow generates a 
complex flow topology due to the highly 3D and unsteady nature of this flow. All 
these features have been resolved in our simulation. Here, we discuss the typical 
features of this flow to demonstrate the accuracy of the current LBE-LES approach. 
Figure 4.3 shows a 3D perspective of a typical flow showing vortical features. It can 
be seen that the horse-shoe vortex at the leading edge of the jet is clearly seen along 
with its interaction with the vortices in the jet shear layer as the jet turns towards 
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COUNTER ROTATING VORTEX PAIR 



Figure 4.3: Instantaneous vorticity magnitude iso-surface |u>| = 0.003. Horse-shoe 
structure (near the wall) and CRVP (following the jet trajectory) are observed (R = 
1 ). 


the outflow direction. Kidney-shaped vortices form along the shear layer, as observed 
earlier [32], Further downstream, the jet is fully oriented towards the outflow and 
contains the classical counter rotating vortex pair (CRVP) seen in nearly all JICF 
experiments. The streamline pattern shown in this figure also shows how the crossflow 
goes around the jet and forms a recirculation region in the wake of the jet. 

At the lee side of jet, the jet deflects by the cross flow and further downstream 
creates a recirculation region. The size of this recirculation zone is about 1.5d for 
R = 1, as measured based on stream-wise velocity at the jet center at z/d = 0.05, 
0.1, and 0.3 and shown in Fig. 4.4. This figure also shows that there is a very small 
reverse flow at the leading edge with the size of 0.25d. This region is related to the 
horse-show vortex that forms on the leading edge of the jet and may be the result of 
an adverse pressure gradient caused by the blockage effect of the jet. 

The origin and formation of the vortices in the wake have been investigated earlier 
experimentally [46]. They examined the direction of wake vortices and horse shoe 
separation at the jet sides. Figure 4.5 shows the normal vorticity at z/d = 0.2 
in the LBE simulation. The wake vorticity at the right side of the jet (shown in 
black) rotates in clock-wise direction and the left side vorticity has counter clockwise 
rotation. The rotational direction of wake vortices is consistent with the experiments 
[46]. The separation of flow at the lee-side of the jet is also clear in this picture. 
Vortex shedding downstream of crossflow is captured as well. 
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Figure 4.4: Normalized streamwise velocity along the stream-wise direction at z/d 
= 0.05, 0.1, 0.3 for R = 1. 



Figure 4.5: Normal vorticity magnitude at the normal cross section plane shows the 
footprints of vorticities at the lee side of the jet (R = 1). 
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(a) Flow through hanging vortex 


(b) Flow in the wake of the jet 


Figure 4.6: Instantaneous streamlines pattern in the side and wake of the JICF. 
Horse shoe streamlines at the lee side of the jet follow a spiral path in the wake of 
the jet (R = 0.5). 


Figure 4.6 shows two sets of instantaneous streamlines. One set originates within 
the jet, and the second set originates from the crossflow boundary layer upstream 
of the jet exit. It can be seen that part of the crossflow is entrained into the wake 
region while some of it is lifted off into the jet shear layer region. There is substantial 
interaction between the jet shear layer and the crossflow boundary layer on the lee 
side of the jet. 

The transient solution of JICF can be used to further highlight the complex flow 
around the CRVP. Analysis is performed for one shedding cycle of the jet with period 
of T = 0.002 (s). This analysis shows how vortices are stretched or realigned in 
different directions and illustrates the origination of CRVP. 

A schematic sketch shown in Fig. 4.7 identifies the vorticity structures in the flow 
field. In this figure, span-wise vorticity at the lee side of the jet and kidney vortices are 
shown by shaded color. In addition, jet side stream-wise vorticity, vertical vorticity 
and CRVP are shown in gray, white and spotted color, respectively. As shown the 
vortex interaction is highly 3D with vorticity in one direction changing into vorticty 
in another direction during the interaction process. 

The steam- wise vorticity is shown in Fig. 4.8 (a), where dark blue indicates the 
positive value (+1) and the light blue marks the negative value (-1). These structures 
appear at the lateral side of the jet due to the gradient of normal velocity (^) at 
the jet boundary layer and the gradient of span- wise velocity (^j) because of the 
higher cross flow deflection near the wall. These structures twist and liftoff from 
the wall region and are torn further downstream (or will change into vorticity in 
another direction). In the past, the location of the vortex breakdown is determined 
by plotting the variation of axial velocity component through a path-line passing 
through a jet side stream- wise vorticity. The location of such breakdown can be 
recognized where there is a “drop” in the velocity component. To identify this, 
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Figure 4.7: Sketch shows the vortex patterns in JICF. 


one streamline (originating from one lateral side of the jet) is tracked through the 
streamwise structure and the axial velocity along this streamline is plotted. 

Figure 4.8 (a) shows the location of this streamline and part (b) shows the axial 
velocity versus distance along the streamline. A sudden “drop” of axial velocity 
appears at distance of 3.5D, which may be used to identify the location of the vortex 
break down. In an earlier study of a circular jet in cross flow [56] this breakdown was 
seen closer to the jet. It is not clear at present if this discrepancy is due to the shape 
of the jet or some other effect. 

The jet side stream-wise structures are lifted into the jet shear layer and therefore, 
are torn by vertical structures along the jet shear layer, as shown qualitatively in Fig. 

4.9. In this figure, jet side stream-wise vorticity is shown in blue and vertical vorticity 
is shown in red. Part (a) of this figure shows how the vertical structures wraps 
around the stream-wise vorticity and in part (b) appear to pinch off the streamwise 
strucure. The detached stream-wise vortex then follows the jet trajectory and creates 
the necessary background circulation for CRVP. 

The vertical structures are also torn during this interaction process, as shown Fig. 

4.10. As can be seen from these figures, there is a complex 3D interaction between 
the shear layers formed along the jet edge and the cross-flow. This interaction results 
in formation of all three components of vorticity and as the jet bends towards the 
outflow direction vorticity in one component changes to another. Vortex stretching, 
breakdown and re-connection of structures occur in this flow field. 
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(a) (b) 

Figure 4.8: Sudden drop in axial velocity indicates the vortex breakdown, (a) jet- 
side stream- wise vorticities, blue: +1 and light blue: -1, (b) Velocity as a function of 
distance along the streamline. 


Other flow structures are also present. Span-wise roller wraps and realigns its 
direction (from span-wise to stream-wise) due to the stronger momentum entrainment 
at the jet side relative to jet center. Figure 4.11 part (a) shows how the surface of the 
span-wise roller is skewed and part (b) shows the resulting CRVP from the realignment 
of span- wise rollers at the corner of the jet. 

In order to confirm the vortex realignment, the contour plots for vortices and 
their stretching components for one shedding cycle of the jet were analyzed and some 
representative figures are shown in Figs. 4.12-4.14. Figure 4.12(a) shows the span- 
wise vorticity at time t = 0.2310 (s) at the jet side ( y/d = —0.5). Also, the span-wise 
vortex stretching component (wjJ^r) and stream- wise vortex stretching component 

(cjj—) are shown in part (b) and (c), respectively. In the region of 0.5 < x/d < 1.0 
and 0.3 < z/d < 0.6, the span-wise vorticity with the magnitude of 1.48 (see part a) 
undergoes a span-wise compression (see part b) and stream-wise stretching (see part 
c). As a result, this span- wise roller at the corner of the jet tilts and combines with 
the torn stream-wise vorticities. 

The span- wise vorticity contours at this time (t + T/2) is shown in Fig. 4.13 (a). 
The tilted vortex has aligned with x-axis and appears as a stream-wise vorticity as 
shown in Figs. 4.13 (b) and (c) in top and side views. In these views, such vortex 
is recognizable at about x/d = 1.0, y/d = -0.5 and z/d = 0.5. The formation of 
the CRVP is seen in the streamwise vorticity at x/D = 1.0. The CRVP is convected 
downstream and its strength changes with time and location. For instance, the CRVP 
at two other locations of x/d = 0.5 and 3.0 are shown in Figs. 4.14 (a) and (b), 
respectively. Signification three dimensionality and interaction of the vortices can be 
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seen in these figures. The lift-off the CRVP is also observable further downstream of 
the jet injection location. 

4.1.4 Mean Flow 

The mean stream- wise velocity profile are presented for R = 0.5 at various stream- 
wise stations ( x/d= 0, 1, and 5) along the jet center plane, y/d — 0 in Figs. 4.15(a-c), 
and along the jet edge plane, y/d = -0.5 in Figs. 4.15(d-f). The current LBE-LES 
results are compared with experimental data and the past numerical predictions using 
various RANS closures [32], At x/d = 0, which is the core of the jet, the crossflow is 
clearly dominant. Further downstream at x/d = 1 a region of recirculation is evident 
in the center plane but is not evident at the edge of the jet, suggesting that the 
recirculation is confined primarily in the shadow of the jet. This is understandable 
since the crossflow is flowing around the jet and provides a confinement effect. These 
observations are consistent with measurements, although the transverse extent (in the 
z-direction) of the reverse flow is somewhat under-predicted by all numerical results. 
Away from the wall (e.g., z/d> 1), the agreement between all numerical results with 
data is quite good. 

Further downstream axially ( x/d = 5), the core flow shows two distinct profiles: 
a near- wall nearly uniform profile for z/d < 1 which may be the jet that has attached 
to the wall and a far field boundary layer like profile for z/d > 1 which may be 
attributed to the crossflow. 

The mean wall-normal velocity profiles along the jet center plane and along the 
jet edge are shown in Figs. 4.16(a-c) and 4.16(d-f), respectively. For this case, the 
jet penetration is limited to z/d < 1 (Fig. 4.16a). Further downstream ( x/d > 1) 
the jet motion normal to the wall is significantly reduced since the jet is now getting 
turned in the axial direction by the crossflow. The agreement in general, is reasonable 
close to the jet orifice and also for z/d > 1 at all locations. There are, however, 
discrepancy further downstream (z/d = 5) in the peak values. The very near- wall 
region ( z/d < 0.25) is not well resolved by all models; however, LBE-LES appears to 
capture the experimentally observed small negative values near z/d = 1 (Fig. 4.16). 
In this region, the flow is turning back towards the wall as the jet is turned down and 
attaches to the wall. 

Figures 4.17(a) and (b) show the computed and experimental contours of the 
normalized mean stream- wise velocity, U/Wj at axial location x/d = 1, respectively 
for R = 1 case. Although the exact shape is not the same the overall features are in 
good agreement, including the small negative velocity in the core. 

4.1.5 RMS and Turbulent Properties 

Turbulent kinetic energy (TKE) profiles along the jet center plane at various axial 
locations are shown in Figs. 4.18(a-c), respectively. In LBE-LES, the TKE is ob- 
tained by time-averaging the resolved kinetic energy and then removing the mean 
kinetic energy, whereas, in the RANS approaches TKE is a modeled quantity. TKE 
is produced not only during interaction between the crossflow boundary layer and 
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Figure 4.12: Vortex tilting at y/d = -0.5, time t. 

(a) Span-wise Vorticity, (b) span-wise stretching, and (c) stream-wise stretching 









- 1.0 - 0.5 0.0 0.5 1.0 


Figure 4.13: Span-wise and stream-wise vorticities at time t + T/2. (a) Span- 
wise vorticity at y/d = -0.5, (b) stream-wise vorticity at z/d = 0.50, (c) stream-wise 
vorticity at x/d = 1.0. ^0 







(b) x/d = 3.0 


Figure 4.14: Stream-wise vorticity at time t + T/2. 






(a) x/d=0 


(b) x/d=l 


(c) x/d=5 



(d) x/d=0 (e) x/d=l 


(f) x/d=5 


Figure 4.15: (a-d) Mean stream-wise velocity profiles (U/Wj) along the jet center 
plane yjd = 0 at x/d = 0, 1, and 5 from the jet center (R = 0.5). (e-h) Mean stream- 
wise velocity profiles ( U/Wj ) along the jet edge plane y/d = —0.5 at x/d = 0, 1, and 
5 from the jet center (R = 0.5). 

— • • RST using Chen (1995) model (Hoda, 2000) 

• Experimental data (Ajersch, 1997) 

— RST using Lam-Bremhorst model (Hoda, 2000) 

— RST using Launder- Tselepidakis model (Hoda, 2000) 

— SRT-LBE-LES using dynamic Smagorisky model 
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(a) x/d=0 (b) x/d=l 


(c) x/d=5 




(d) x/d=0 (e) x/d=l (f) x/d=5 

Figure 4.16: (a-c) Mean wall-normal velocity profiles ( W/Wj ) along the jet center 
plane y/d = 0 at xjd = 0, 1, and 5 from the jet center (R = 0.5). (d-f) Mean wall- 
normal velocity profiles ( W/Wj ) along the jet edge plane y/d = —0.5 at xjd = 0, 1, 
and 5 from the jet center (R = 0.5). 
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(a) LBE Results (b) Exp. data 

Figure 4.17: Comparison of LBE results with experiment for R=l. Normalized 
stream- wise velocity at x/d = 1. 



the jet, but also due to various velocity gradients and strong streamline curvature. 
Peak of TKE profile is located at the location of maximum velocity gradients. TKE 
is larger at the top of the jet where the interactions are dominant (see discussion 
on coherent structures above). Comparison (not shown) of turbulence fluctuation 
intensity in the stream-wise (u rms /Wj) and wall-normal (w rms /Wj) directions shows 
that these fluctuation levels are similar in magnitude but that there are significant 
anisotropy in the near- wall region. This is to be expected since in this region vortical 
breakdown occurs in a complex 3D manner. 

Figures 4.19 (a-c) shows the normalized shear stress (liiv/Wj) profiles along the 
centerline at various axial locations. Although the general trend is similar to the 
experiment, the LES-LBE prediction is significantly larger in many regions. There 
are two possible reasons for this. In general, the Reynolds stress prediction in highly 
separated flow regions is very difficult and the algebraic subgrid closure used here may 
not be appropriate. Another possible reason is that in LES higher order moments 
typically take much longer to reach statistically stationary state and the current 
simulation may not have been run long enough (although around 4 flow-through- 
times have been averaged for the current data). Improved non-equilibrium subgrid 
model, such as the localized subgrid kinetic energy model [24] and longer statistical 
averaging may improve this prediction and will be addressed in the future. 

The comparison between LBE-LES and experiment of the cross-stream (y — z) 
contours of the mean TKE and Reynolds stress (uw/W'j) are shown in Figs. 4.20(a,b) 
and Figs. 4.20 (c,d), respectively. For direct comparison, the same contour interval is 


34 





(a) x/d=0 (b) x/d=l (c) x/d=5 

Figure 4.18: Mean turbulence kinetic energy y/k/Wj along the jet center plane 
y/d = 0 at x/d = 0, 1, 5 from the jet center (R = 0.5). 



(a) x/d=0 


(b) x/d=l 


(c) x/d=5 


Figure 4.19: Turbulence shear stress uw/Wf profiles along the jet edge plane y/d = 
0 at x/d = 0, 1, and 5 from the jet center (R = 0.5). 
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used. In general, there is reasonable agreement between the two studies. The location 
and magnitude of both TKE and Re-stress are well predicted in the current LBE-LES 
study. The normalized contours of the Re-stress vw/Wj shown in Figs. 4.20 (e,f) 
compares the predictions with measurements. Again, there is good agreement with 
data. 

4.1.6 Comparison of SRT and MRT Results 

To avoid oscillations near the sharp corners due to the geometrical singularities at 
these locations, the multiple relaxation time method has been developed in the past 
[32]. Such method eliminates the strong local spatial oscillations and improves the 
stability of the LBE scheme. In the past, comparison between SRT and MRT results 
are done for Stokes first problem, steady uniform flow over finite/semi-finite plate and 
cavity flow [32]. These results suggest that MRT model is better behaved in flows 
involving large gradients and high Re numbers. The Re number for the jet-in-cross 
flow simulation presented here, is about 4700 and large oscillations were not seen 
for SRT simulation. In fact, as shown in the previous section, the SRT model is 
quite capable of capturing most of the features of interest in the JICF for the cases 
simulated here. However, in order to ensure the accuracy of the results, two cases of 
R = 0.5 and 1.0 were also simulated using MRT model and results were compared 
with SRT as well as experimental values. No major difference between SRT and MRT 
results were seen. Some results and included here for completeness. 

Results at x/d = 1 along the jet center and jet edge planes are shown in Fig. 4.21. 
This location ( x/d = 1) was chosen for comparison since it was at the vicinity of the 
sharp edges and the most likely place where results could differ. Figures 4.21(a) and 
(b) show respectively, the stream-wise velocity at the jet center and edge planes, and 
compares them to the SRT results. The scatter between MRT and SRT is negligible 
for this test case. 

Figure 4.21(c) and (d) show the wall normal velocity at center and edge jet planes 
and Figs. 4.21(e) and (f) show the mean turbulent kinetic energy and uw shear stress 
along the jet center plane. Again, it can be seen that there are not that significant 
differences seen between the two methods, at least for this test case. 

Nevertheless, based on past studies, it is very likely that the MRT LBE model is 
probably the more robust one and should be used for high Re test cases. 

4.2 Hybrid FV-LBE-LES Solver Validation 

In this section, we first test the boundary information exchange in the coupling be- 
tween the two solvers by passing an acoustic pulse from each domain, and then de- 
scribe the coupled simulation of JICF case discussed above to complete the validation 
process. 
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(c) LBE Results (d) Exp. data (e) LBE Results (f) Exp. data 

Figure 4.20: Comparison of LBE results with experiment for R = 1. (a,b) Nor- 
malized TKE ^ at x/d — 8, (c,d) Normalized shear stress at x/d = 5, (e,f) 
Normalized shear stress ^ at x/d = 5 
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(a) U/Wj (b) U/Wj (c) W/Wj 





(d) W/Wj (e) Vk/Wj (f) uw/Wj 


Figure 4.21: Multiple-relaxation-time results for R = 0.5 shown with triangles on 
top of single-relaxation time results shown with solid lines. (a,b) Mean stream-wise 
velocity along the jet center and edge planes, (c,d) Mean wall normal velocity along 
the jet center and edge planes, (e) Mean turbulent kinetic energy along the jet center 
plane, (f) Turbulence shear stress profiles along the jet edge plane. 
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4.2.1 Acoustic Pulse Test 


The behavior of the coupling scheme at the boundaries is examined using acoustic 
wave test. This type of test is used in many codes to evaluate the proper implemen- 
tation of boundary conditions. In the first case of interest, the LBE to FV coupling is 
tested. An acoustic pulse, is initially generated in LBE domain, propagating toward 
FV domain. The acoustic pulse was is created at t = 0 by using the initial conditions: 


u- = u; + Aexp (b ( V* L f /2 ) 

<?(!)=! + W-U-) 
where Uq = 0, A = 0.01, and B = 15 are considered. Figure 4.22 shows the velocity, 
density and pressure variation along the stream-wise direction. The amplitude and the 
wave length of the pulse remains constant as it passes through the coupling boundary 
which proves the validity of the coupling approach when information are sent from 
LBE to NS domain. 

The characteristic velocity for LBE is 0.05 and is equivalent to 10.15 m/sec for FV 
formulation. Also the characteristic length for LBE is 80 which is equivalent to 0.01 
(m) . Equivalent velocity and length for FV domain are determined by matching the 
Re number. The computational grid for LBE and FV domains are 80X20X20 each. 
This case is setup in such a way that time step for the LBE is equal to time step for 
FV scheme. The time step for LBE scheme is computed as: 



Cr 


0.000125 _ 
203 


5t les = CFL *{) = 0.5 x 10" 6 

Finally, the transfer from the FV to the LBE is tested by generating the acoustic 
wave inside of the FV domain, which propagates toward LBE domain. The initial 
conditions at t = 0 are set as: 


u = uq + A exp I — ( B 


, a:-L/2 y 


P = Po + PqCq(u — Uq ) 

„ „ . Po(u-u 0 ) 

P = P 0 + 


Co 


T = 


pR , 


: gas 


where u 0 = 10.15 (m/s), A = 4, B = 6, p 0 = 101325(po), Co = y ^yR^To, 7 = 1-4, 
T 0 = 300(A), R gas is the gas constant. 

Figure 4.23 shows the velocity, density and pressure variations along the stream- 
wise direction. The amplitude and the wave length of the pulse remains constant as 
it passes through the coupling boundary which proves the validity of the coupling 
approach when information is sent from the FV to the LBE domain. 
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LBE Part LES Part 

Figure 4.22: Acoustic pulse from the LBE domain to the FV domain 
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LBE Part LES Part 

Figure 4.23: Acoustic pulse from the LBE domain to the FV domain 
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4.2.2 Single Jet in Crossflow 

In this study, the experiment of Ajersch et al. [51] is again chosen as the benchmark 
case for validation. The dimensions of the computational domain are shown in Fig. 
4.1. The simulation is carried out at Reynolds number of 4700 (based on the jet 
velocity and the nozzle width D ) and at a jet-cross-flow velocity ratio of 0.5. The 
cross-flow velocity profile is initialized with a boundary layer thickness of 5 D. The 
computational domain is resolved using 128 x 96 x 64 for the cross-flow domain (FV- 
LES) and 20 x 20 x 100 for jet section (LBE). So a total of 0.82 x 10 6 grid points 
is used to discretize the complete domain. The results from these simulations are 
compared with the earlier results obtained using the LBE-LES for the entire domain 
[ 12 ]. 

Computation of a full simulation requires approximately 1155 single-processor 
hours with 1.5GB memory on the IBM SP4 machine. Same boundary conditions as 
for single JICF in LBE simulation were specified here as well. 

Results are compared with the experiment of Ajersch and the earlier full LBE-LES 
simulation. Comparison of LBE-LES with numerical results of Hoda et al [52] were 
discussed earlier and it was shown that the full LBE-LES approach can capture with 
reasonable accuracy the experimental mean and rms velocity profiles. Therefore, we 
only compare the current LBE-FV-LES coupled results against the earlier LBE-LES 
results. 

All the important features in JICF such as the horse-shoe (or kidney-shaped) 
structure and counter rotating vortex pair (CRVP) observed in earlier LBE-LES (and 
in the experiments) are captured in the present simulation with reasonable accuracy 
(and in very good agreement with the previous study). 
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Figure 4.25: Velocity contours shows the flow in the recirculation zone and the 
motion downstream. 


Mass conservation of the inflow and outflow is shown in Fig. 4.24 to demonstrate 
that this coupled simulation approach can conserve mass. The velocity contours in 
Fig. 4.25 shows the formation of the recirculation zone behind the jet. By reviewing 
time evolution of the flow field it is seen that the recirculation bubble periodically 
sheds the flow in-phase with the motion of the kidney-shaped vortices observed in 
this flow. 

The periodic shedding of the flow from the recirculation zone is show more clearly 
in Fig. 4.25 which shows the velocity contours along stream-wise direction. Smooth 
transition of velocity contours from LBE domain to NS domain shown in this figure 
clearly demonstrated the accuracy of the LBE-NS coupling employed in the present 
approach. 

The mean velocity profiles comparison with past LBE-LES and experimental data 
at representative stream- wise stations ( x/D = 0,1,3) along the jet center plane 
(Y/D = 0) is shown in Figs. 4.26 (a-i). The exit plane velocity is well captured 
by LBE-LES and the reverse flow at X/D = 1 is also captured. Turbulent profiles 
at the jet center plane are shown in Figs. 4.27(a-i) which also shows good agreement 
with experiment. 

4.2.3 Parallel Performance and scalability of the Hybrid Solver 

The overall performance of the coupling algorithm is determined for the above JICF 
case. The test grid resolution of 72x64x60 & 25x25x50 for crossflow and jet sections 
are considered in this study. So the total of 300k grid points, which require 1.07 
Gbytes memory are solved on a IBM/SP4 platform to obtain this scaling data. 

Table 4.2 lists the timing results for four cases with different number of processors. 
In this table, the total wall time per step, scalability and the ratio of number of grids 
in each processors of FV domain to the number of grids in each processors of LBE 
domain are given. 
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(a) (U/Wj) at (b) (U/Wj) (c) (U/Wj) 

x/D—O at x/D— 1 at x/D= 3 


X / D = 3 



(d) {V/W 5 ) at (e) (V/Wj) (f) V/Wj) at 

x/D—O at x/D= 1 x/D= 3 



(g) TO) at (h) (W/Wj) (i) (W/W, ) 

x/D=0 at x/D= 1 at x/D= 3 


Figure 4.26: Mean velocity profiles along the jet center plane (y/D= 0) at x/D= 0, 
1, and 3 from the jet center. 
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(a) 

(^rms /Wi) (b) (c) 

at x/D = 0 (fVms/ttj) {Uj •ms/Wj) 

at x/D = 1 at x/D = 3 



(d) (w rms /Wj) (e) (f) 

at x/D = 0 {Wrms/Wj) (Wrms/Wj) 

at x/D = 1 at x/D = 3 



(g) ( uw/Wf ) (h) (i uw/Wf ) (i) (i uw/W ? ) 

at x/D = 0 at x/D =1 at x/D - 3 


Figure 4.27: Turbulent profiles along the jet center plane (y/D= 0) at xj D= 0, 1, 
and 3 from the jet center. 
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Case # 

1 

2 

3 

4 

Number of CPU’s (FV + LBE) 

5 (4+1) 

9 (8+1) 

17 (16+1) 

31 (30+1) 

Run-time/step (secs) 

0.76 

0.38 

0.25 

0.23 

Scalability (T ba3e /T N ) 

1.0 

2.0 

3.0 

3.3 

Grid ratio (FV/LBE) 

3.0 

1.5 

0.7 

0.4 


The scalability is defined as T base /T N where T base is the run time on N base processors 

and T n is the run time on N processors. 


Table 4.2: Coupled scheme performance and scalability on a IBM/SP4 for the single 
JICF validation case with a resolution of 300K grid points. Nha Se —4 and T base = 12.8 
(min.). 


Case # 

3 

Number of CPU’s 

17 

Run-time/step (secs) 

1.21 

Parallel cost(% of run time) 

14.9 

Cost for FV domain (% of run time) 

80.1 

Cost for LBE domain (% of run time) 

5.0 

Scalability 

3.0 


Table 4.3: Computational cost for the single JICF validation case with a resolution 
of 300K grid points. 


A performance analysis for case #3 is carried out in order to measure the parallel 
overhead and the cost of each subroutine in the coupled algorithm. Table 4.3 summa- 
rizes the percentage of run time for the Navier-Stokes solver, LBE solver and parallel 
communication. As expected FV domain takes the majority of the run time, i.e. 80%, 
since it solves 7 finite-volume differenced equations on a large grid. The LBE only 
solves for the advection equation (along 18 directions) on a relatively small grid size 
with a cost of 5 percent of run time. The parallel cost is large at about 15% due 
to the communication time between the two domains, as well as the communication 
time for the solvers themselves. 


4.3 Multiple Micro- Jets in Cross-flow 

This section reports on the fully coupled LBE-FV-LES simulation of turbulent mul- 
tiple micro-jets in crossflow. Earlier micro-hole matrix of different combinations were 
simulated. Here, we report on a representative case of a matrix of 3 x 3 micro holes 
in a high Mach number turbulent boundary layer. This test case is similar to the 
experiment at NASA/GRC [4] except that only a 3 x 3 matrix is simulated, whereas 
in the experiment, thousands of micro holes on a large flat plate were employed. Our 
present effort is to demonstrate this coupled approach. 

In the micro blowing technique experimentally investigated at NASA/GRC [4], a 
plate with porosity of 13 to 43 percent was studied and it was shown that drag can be 
significantly reduced. Holes diameters ranged from around 0.2 to 0.5 millimeter and 
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Figure 4.28: Geometry and grid configuration for multi hole simulation (every five 
grid is shown). 


depending on plate porosity, an order of 100000 micro holes existed on a 100cm 2 flat 
plate area. Simulation of flow inside such tiny holes along with a cross-flow domain 
whose length scale is orders of magnitude larger is a challenging problem. 

Here, we simulate a 3 x 3 square pattern of injection holes each with diameter 
of D = 0.5 mm and equally spaced a distance of D apart, as shown in Fig. 4.28. 
Turbulent inflow profile is specified 5 D upstream of injection holes with a free stream 
velocity of Uo = 138.54 m/s and boundary layer thickness of 5 = 12 mm. These 
conditions closely approximate the experimental data. Slip boundary conditions is 
imposed on the upper plane at a distance 6/2 above the injection holes, and charac- 
teristic outflow conditions were imposed at a location 12 D downstream of the holes. 
Periodic boundary conditions are imposed in the spanwise direction. 

The plate porosity is defined by cumulative hole area divided by the plate area. 
For our test case here it is 25%. The blowing rate (Vj et /Uo)) of 0.006 is simulated 
here. The holes have an aspect ratio ( L/ D ) of 4 and the hole inflow velocity profile 
is specified as: 

Tr T , . ,7TX. . ,7 TV. 

Vinflow — ) • 

Based on the experiments done at NASA/GRC, the flow conditions are set as 
following. The tunnel pressure is 0.24 atm and the density, the kinematic viscosity 
and the Reynolds number per meter are found as, 


P = 


P 

rT 


24318 
(287) (300) 


0.28 (kg/m 3 ), 


a 18.48 x 10~ 6 „„ , . , 

" = ~p = 0^8 = 66 '° X 10 " (m /s) 


Re = 


Vx Re_ 138.0 
v x 66.0 x 10 6 


= 2.1 x 10 6 (/m) 
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Figure 4.29: Inflow boundary layer profile imposed for numerical simulation which 
matches the inflow BL measured in the experiment. 


Since the leading edge of the test plate in the experiment is 25.4 cm downstream from 
the end of wind tunnel transition duct, the Re number at the leading edge of the test 
plate is found as, 

Xstart — 0.254(m) =>■ Re x = 0.55 x 10 6 
Using this information, the momentum thickness Re can be calculated, 

Re e = OM5Re x 079 = 2683 


which gives the momentum thickness of 9 = 0.0012m. The boundary layer shape 
factor measured in the experiment is equal to H = 1.34, and thus the theoretical 
friction coefficient is computed [57] as, 


C f = 0 . 246. 1 0 _0 - 678/f ( ) ~ 0 - 268 = 0.0039 

H 


In addition, the theoretical boundary layer thickness can be determined from the 
shape factor and the momentum thickness using the following relation, 

H = l + - = 1.34 => n = 5.88 

n 


e 

8 


n 

( n + 1 )(n + 2) 


0.1084 


5 = 0.0126 


Using this information, a turbulent inflow boundary layer is created for the nu- 
merical simulation and shown in Fig. 4.29. In this graph the theoretical boundary 
layer thickness is highlighted. 
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Figure 4.30: Time averaged velocity profile at the inflow based on the wall coordi- 
nate. 


This inflow is generated in such a way that agrees with the Log-law profile as well. 
Figure 4.30 depicts the time averaged velocity profile based on wall coordinate. The 
log law profile is shown in this graph as well for comparison. 

Figure 4.31 shows the boundary layer profile at different stream-wise locations 
along the plate in-between the injection holes. It can be seen that even this small 
blowing can result in the change in the boundary layer profile in the inner layer 
(y/S << 0.03). 

Since the injection starts at x/D = 5, the friction velocity decreases past this 
location, and once the injection stops (at x/D = 11), it begins to recover to its 
original value. This is in good agreement with experimental observation. Figure 4.30 
shows the comparison of the experiment profile at the inflow with the current inflow 
profile. It also shows the velocity profile at the injection port. As we can see, since 
the friction velocity goes down due to injection, higher U + is achieved, as expected. 

Velocity vector field in the center plane for the three injection holes are shown in 
Fig. 4.32. In the microblowing approach, the flow from the injectors is very small 
and with a very small velocity. Thus, these jets interact with the viscous sublayer 
and the buffer layer of the incoming crossflow boundary layer. Analysis shows that 
the interaction is limited to a region up to y + = 19. It can be seen from this figure 
that whereas the jet from the leading injector penetrates further into the crossflow, 
the jets in its shadow have weaker impact and this shadow effect increases for the 
third injector. 

Further analysis shows that there is a complex interaction between the various 
injectors. Results suggest that the wake effect of the leading injector can change the 
local pressure above the injectors behind it and can cause periodic reduction of the 
mass flow. More detailed analysis is still needed to fully understand the dynamics of 


48 





Figure 4.32: Velocity field in the center place of the three injectors. The shadow 
effect of the leading injector is seen on the injectors behind it. 


this interaction process as a function of the injection pressure. 

3D visualization of the coherent structures in the near-wall region is shown in Fig. 
4.33 in the y + = 4 — 19 region. The flow is highly modulated by the presence of 
the micro injection and only further downstream do we see the formation of more 
turbulent mixed flow. Streaks in the near wall region can also be seen downstream 
of the martix of the holes. 

Spanwise-averaged normalized coefficient of friction in stream-wise direction for 
two different injection ratio of 0.02 and 0.07 and shown in Fig.4.34. As observed in 
the experiments [4] we also observe a significant drop in skin friction drag by up to 
50 percent in the immediate vicinity of the injectors. 

Finally, Figs. 5.1 shows some of the flow features near the injector holes showing 
that there is significant interaction between the injectors. Interaction increases as the 
flow moves away from the plate and most of the interaction occurs in the viscous sub 
layer and in the buffer layer. Since most of the TKE production is in the buffer layer, 
this interaction results in a major modification of the near wall boundary layer profile 
and hence, impacts the wall shear stress and the resulting skin friction drag. 


L 
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CHAPTER V 


CONCLUSIONS AND FUTURE PLANS 


This final report describes the development and the evaluation of the fully coupled 
LBE-LES and FV-LES simulation of turbulent jet in crossflow and multiple micro-jets 
in crossflow. Earlier, the single free jets and jet in crossflow (JICF) were simulated 
using a LBE-LES method [12, 13, 10, 14, 11] and compared to experimental data. 
Very good agreement was obtained. Here, we simulate the same test case using the 
coupled formulation and compare to earlier studies to demonstrate the applicability 
of this coupled method. Subsequently, the coupled solver is used to simulate a matrix 
of 3 x 3 micro holes in a high Mach number turbulent boundary layer. This test case 
is similar to the experiment at NASA/GRC [4, 6] except that only a 3 x 3 matrix is 
simulated, whereas in the experiment, thousands of micro holes on a large flat plate 
were employed. 

This is the first reported development of a coupled LBE-FV-LES methodology for 
application to micro-blowing technique. The coupled LBE-FV-LES can also be used 
for other possible scenarios such as film cooling injection, distributed control by fuel 
injection, etc. 

It has been shown that this coupled approach can handle many holes simultane- 
ously in a cost effective manner. We also believe that these sort of sub-set simulations 
can be used to investigate similarity features and also to develop new scaling laws. 
However, due to resource and time constraints, some of the planned tasks had to be 
truncated. Nevertheless, as this report demonstrates, we now have a hybrid LBE-FV 
capability to carry out LES of this complex problem. 

In the current study, the LBE-LES is employed to simulate the flow inside the jet 
nozzles while the FV-LES is used to simulate the crossflow. A single jet in crossflow 
case is used for validation purpose and the results are compared with experimental 
data and full LBE-LES simulation. Good agreement with data is obtained. Subse- 
quently, MBT over a flat plate with porosity of 25% is simulated using 9 jets in a 
compressible cross flow at a Mach number of 0.4. It is shown that MBT suppresses 
the near-wall vortices and reduces the skin friction by up to 50 percent. This is in 
good agreement with experimental data. 

Currently, there are no immediate plans to further study this MBT problem since 
funding for this program has ended. However, it is planned to keep this code active 
and further improved for performance. Additional changes being considered for im- 
plementation in this code are (a) an ability to simulate circular or elliptical holes, (b) 
an allowance for an angled injection and (c) a more generalized subgrid closure for 
the LBE using the LDKM model. Some progress has already been made to achieve 
these objectives but more effort in validation and testing is still needed. 
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(e) y+=6.82 (f) y+=8.09 

Figure 5.1: Normal velocity at different wall normal locations for R = 0.003 
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